Efficient generation of correlated random numbers 
using Chebyshev-optimal magnitude-only IIR filters 
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We compare several methods for the efficient generation of correlated random sequences (colored 
noise) by filtering white noise to achieve a desired correlation spectrum. We argue that a class of 
IIR filter-design techniques developed in the 1970s, which obtain the global Chebyshev-optimum 
minimum-phase filter with a desired magnitude and arbitrary phase, are uniquely suited for this 
problem but have seldom been used. The short filters that result from such techniques are crucial 
for applications of colored noise in physical simulations involving random processes, for which many 
long random sequences must be generated and computational time and memory are at a premium. 
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The generation of correlated random sequences, or 
"colored noise" , is important for many physical simula- 
tions involving random processes 0, 0,11, 0, 0] , and often 
the required computational time and memory is a criti- 
cal concern. Although many filter-based technicmes have 
been applied to this problem [l|, 0, 0, 0, 0, 0, El , in 
this brief manuscript we point out that these methods are 
suboptimal: a global Chebyshev-optimum stable IIR fil- 
ter for this problem may be designed based on techniques 
developed in the 1970s [13, El- 

Colored-noise generation is required for many types of 
numerical simulations, such as thermodynamics 0,0,111], 
laser noise and first-passage time problems [if, and 
chaotic dynamics [Bj]. In general, any numerical model 
involving stochastic differential equations in which there 
is some background distribution, nonlinearity, or external 
noise associated with the quantities driving the fluctua- 
tions will require the use of colored noise 0, 0, 0, 0, E3] • 
Such colored noise is most commonly described by a lin- 
ear process: the output of a time- invariant linear filter 
applied to white noise 14[. The key implementation 
question, as we discuss below, is to determine how to 
apply this filter process. With rare exceptions [H|], the 
underlying random process is almost always Gaussian in 
nature, if it arises by a physical process governed by the 
central-limit theorem |16j |. It is this nearly ubiquitous 
case that we focus on here, in which only the correlation 
function of zero-mean colored noise is of concern, and not 
any higher-order statistics. 

From a computational standpoint, the central problem 
is that many of these applications require many random 
sequences to be generated in parallel, and/or very long se- 
quences, imposing performance and memory constraints 
on the noise generation. This occurs, for example, in 
the generation of power-law (1/ f a ) "pink noise" for ap- 
plications ranging from hydrology to music, where very 
long sequences (whose length may not be known a priori) 
are often required Multiple long random sequences 
also occur in the generation of Rayleigh random vari- 
ates for simulation of nonisotropic scattering, macrocel- 
lular systems, and physical models of mobile multipath 




FIG. 1: Plot of Planck distribution R(uj) — au> / (exp(aa;) — 1) 
vs. lu (red). The periodogram (spectrum) of a finite-length 
correlated random sequence generated by the methods in 
this paper is also shown (black), computed by Bartlett's 
method [18f |. 



radio channels [3, U, 17 1. Perhaps the most stringent re- 
quirements for colored-noise generation occur for spatio- 
temporal noise: where a different random sequence may 
need to be generated simultaneously for every point in 
space. This occurs, for example, in modeling of Brownian 
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motion [B[ , or for simulation of thermal radiation 

The fundamental technique for generating correlated 
random sequences is to start with white noise (uncor- 
rected random numbers) and to apply a filter whose 
power spectrum matches the desired correlation spec- 
trum [3, 0, 0, EU 11 1- That is, suppose that we want 
colored-noise y n with some desired correlation function 
Rm = DnDn+m, the (discrete-time) Fourier transform of 
which is the correlation spectrum R(cu). We then start 
with white noise x n , usually uncorrelated Gaussian ran- 



dom numbers, with zero mean and 



1, and apply 
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a filter H(uo): in frequency domain, Y(u>) = H(u>)X(uj). 
The desired correlation spectrum is achieved if \H (w)| 2 = 

Although this filtering operation can be performed en- 
tirely in the frequency domain via a fast Fourier trans- 
form (FFT) of the data sequences 0, 0, such an 
approach is often too inefficient. Filtering entirely in the 
frequency domain requires the entire data sequence y n 
to be computed and stored in advance, and if many long 
sequences are required the storage becomes prohibitive. 
The alternative is to perform the filtering in the time do- 
main, using either finite-impulse response (FIR) filters or 
infinite-impulse response (IIR) filters (also called recur- 
sive or ARM A filters). Since the generation of colored 
noise is of interest to many researchers without a back- 
ground in signal processing, let us briefly review the ba- 
sics of FIR and IIR filtering, which are covered in more 
detail by numerous textbooks, e.g. Ref. [H- FIR filters 
consist of a finite-length sequence b n that is convolved 
with the x n . The Fourier transform of the b n is the fil- 
ter H(lu), but in general a finite-length sequence can only 
approximate an arbitrary desired spectrum R(uj). In par- 
ticular, an FIR filter yields a spectrum H(u>) which is a 
polynomial in z = e luJ with coefficients b n . A better 
approximation may be obtained by generalizing to a ra- 
tio of two polynomials (i.e., rational functions), which 
leads to IIR filters. An IIR filter is determined by two 
finite-length sequences, a n (n = . . . N, ao = I) and 
b n (n = . . . M) , that determine y n via the recurrence 
(which can be written in several equivalent forms): 



optimization 11 1 or Yule- Walker methods [6[ that are 
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That is, y n is a convolution of b n with x n and of a n with 
the previous values of y n . The filter design problem is 
then, given filter orders N and M, to find the a n and b n 
that best approximate the desired spectrum. 

Therefore, we must choose what type of filter to ap- 
ply (IIR or FIR) and a filter-design method. Several 
choices have been previously proposed in the context of 
colored-noise generation. The simplest method, as we 
mentioned above, is to just perform a fast Fourier trans- 
form (FFT) of the entire desired spectrum multiplied by 
random phases 0, H, @, [l(|. This is equivalent to an 
FIR filter of the same leng th as the data, designed by 
the "window" method 18]. One can also employ FIR 
filters of shorter lengths, designed by a variety of stan- 
dard methods such as Parks-McClellan [3, [l^] . For cer- 
tain noise problems, an FIR filter may also be designed 
by analytical methods [1, 0]. In order to shorten the 
length of the required filter, and thus the memory and 
time requirements, IIR recursive filters have been pro- 
posed [H, E| EH- I n general, the design of IIR filters 
is a difficult problem [18|, |20(, and past approaches to 
colored noise generation by IIR filtering have used local- 



not guaranteed to yield the global-optimum filter coeffi- 
cients. Another difficulty is that not all IIR-filter design 
techniques are guaranteed to yield a stable filter (one 
which does not lead to a diverging process). One im- 
portant exception is exponentially correlated noise, for 
which an exact first-order IIR filter is known analytically 
and is commonly used (although typically derived from 
a stochastic differential equation y' = —ay + x and not 
recognized as an IIR filter per se) [J SHU- 

However, 

there is a key property of the IIR filter design problem for 
colored-noise generation that makes optimal filter design 
practical: the phase of the filter H (cu) is irrelevant, since 
it is multiplied in any case by white noise X(lj) with a 
random phase. 

In particular, we can exploit results by Dud- 



geon [12[ and Rabiner [13|, who demonstrated that the 



global Chebyshev-optimum magnitude-only minimum- 
phase stable IIR filter-design problem can be effi- 
ciently solved by a sequence of linear-programming prob- 
lems 2l|. Alternative methods with similar properties 
have also been proposed [22I , |23| , and the Dudgeon and 
Rabiner technique was recently generalized to multidi- 
mensional IIR filters 2J|. However, its applicability to 



the problem of colored-noise generation does not seem to 
have been appreciated, and in this manuscript we demon- 
strate that it can yield dramatically more efficient filters 
than previous approaches. 

To demonstrate the efficacy of the various filter-design 
approaches for colored-noise generation, we consider an 
example drawn from thermodynamic simulations of gray- 
body thermal radiation 0]. In this problem, thermal ef- 
fects are modeled as random fluctuating current sources 
everywhere in space, with a correlation spectrum R(u) = 
auj / (exp(auj) — I), for a constant a determined by the 
temperature, based on the Planck distribution. This dis- 
tribution is shown in Fig. I, along with the periodogram 
of a finite-length correlated random sequence generated 
by the methods in this paper. 

An IIR filter is defined by a rational polynomial func- 
tion in z = e tu : 
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where N and M define the filter order and ao — 1 for 
convenience. A stable minimum-phase IIR filter has all 
of its poles and zeros within the unit circle in z (i.e., for 
Imw > 0) [18]. Given a sequence of uncorrelated random 
numbers x n , the output colored-noise sequence y n is then 
given by the recurrence (Eq. 1) above. The required 
memory, along with the computation time per output, 
is therefore 0(N + M). An FIR filter is the special case 
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FIG. 2: Plot of Chebychev error, = max^ \R(uj) — 

\H(u>)\ 2 \, vs. N + M for the Planck distribution R(uj) 
given in Fig. 1. The Chebychev error is plotted for six dif- 
ferent filter methods: IIR Yule- Walker (red squares), FIR 
Parks- McClellan (black crosses), FIR Least-Squares (blue 
dots), IIR global-optimum filter (gray circles), and two non- 
linear conjugate-gradient methods optimizing two different 
L 2 norms: (a) L2(|//| — VR) (green diamonds) and (b) 
L2(\H\ 2 — R) (magenta triangles). 
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FIG. 3: Plot of L 2 error, L 2 {\H\ - VR), vs. N + M for the 
Planck distribution R(u) given in Fig. 1. The L 2 (\H\ - y/R) 
error is plotted for six different filter methods: IIR Yule- 
Walker (red squares), FIR Parks- McClellan (black crosses), 
FIR Least-Squares (blue dots), IIR global-optimum filter 
(gray circles), and two nonlinear conjugate-gradient methods 
optimizing two different L 2 norms: (a) L 2 (\H\ — VR) (green 
diamonds) and (b) L 2 (\H\ 2 — R) (magenta triangles). 



N = 0. A sequence y n of length K from an FIR filter can 
be generated in 0(K log M) time instead of O(KM) by 
use of fast Fourier transforms (via overlap-add or overlap- 
save techniques fl8j|). but the memory requirements are 
not improved. 

In Fig. 2, we plot the (Chebyshev) error, 

max w \R(lj) — \H(uj)\ 2 \, as a function of the total fil- 
ter order N + M, for filters H(uS) designed by several 
techniques. (Here, we employ the Chebyshev error over 
the entire frequency bandwidth; for other applications, 
only a subset of the bandwidth may be of interest.) 
The best method, i.e. smallest Chebychev error for 
any given-order, is the optimal magnitude-only IIR fil- 
ter design (implemented using the differential-correction 
algorithm [l2l . l2l]|). The other design methods plot- 
ted consist of two FIR and two IIR filter techniques. 
The two linear-phase FIR filter techniques are (from 
the Matlab signal-processing toolbox): first, the Parks- 
McClcllan algorithm, which finds the global Chebyshev 
optimum fT8l. Tl9l] (black); second, a least-squares FIR op- 
timization (blue), as described in Ref. 0. The two IIR 
filter designs plotted are: first, a nonlinear conjugate- 
gradient minimization of a least-s quar es norm proposed 
by Ref. HH, and suggested by Ref.lUI for use in genera- 
tion of Rayleigh-correlated noise (green diamonds); sec- 
ond, another (non-global) optimization technique based 
on the modified Yule- Walker algorithm in the Matlab 



signal- processing toolbox [6|, [26[. The conjugate-gradient 
method only finds a local optimum, and is highly sen- 
sitive to the starting point of the optimization for this 
problem Here, we use the Yule- Walker IIR filter as 
the starting point, and conjugate gradient is able to im- 
prove upon it significantly (in fact, it happens to nearly 
match the Parks-McClellan FIR performance). We ex- 
amined the conjugate-gradient minimization of two dif- 
ferent I/ 2 norms [L 2 (x n =i... n) = \/J2 n X n/ N \ °f tne er_ 
ror: L 2 (\H\ - y/R) 0, SI and L 2 {\H\ 2 - R), although 
it turns out to make little difference in the result. 

We suspect that the Chebyshev norm is typically the 
most appropriate one for physical simulations involving 
random processes. The reason is that a large error in a 
narrow bandwidth, which might be allowed by a least- 
square (L 2 ) norm, could result in a spectral feature that 
might be mistaken for a spurious physical phenomenon; 
a large "spike" in error may also interact adversely with 
nonlincarities in the physics. Nevertheless, in Fig. 3 and 
Fig. 4, we show two different L 2 errors, L 2 (\H\-y/R) and 
L 2 (\H | 2 — R), for the same filter designs as in Fig. 2, and 
the results demonstrate that the Chebychev-optimal IIR 
filter is at least comparable, and often superior to, the 
other methods, even in norms that it does not strictly 
optimize. In particular, the L 2 (\H\ 2 — R) norm of the 
Chebyshev-optimum IIR filter does better than all other 
local optimization techniques for most N + M. Further 
gains could potentially be made, if this is the desired error 
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FIG. 4: Plot of L 2 error, L 2 (\H\ 2 - R), vs. N + M for the 
Planck distribution R(u) given in Fig. 1. The L 2 (\H\ 2 - R) 
error is plotted for six different filter methods: IIR Yule- 
Walker (red squares), FIR Parks- McClellan (black crosses), 
FIR Least-Square (blue dots), IIR global-optimum filter (gray 
circles), and two nonlinear conjugate-gradient methods op- 
timizing two different L 2 norms: (a) L 2 (\H\ — \fR) (green 
diamonds) and (b) L 2 (\H\ 2 — R) (magenta triangles). 



norm, by starting with the Chebyshev filter and then 
performing a local optimization by some other method. 
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